home *** CD-ROM | disk | FTP | other *** search
/ Deutsche Edition 2 / Deutsche Edition 2.iso / mac / UTILITYS / MathLibFast ƒ / MathlibFast.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-14  |  3.4 KB  |  180 lines  |  [TEXT/MPS ]
  1. /*
  2.     source code for replacement for math library
  3. */
  4.  
  5. #include <math.h>
  6.  
  7. #define    pi_ 3.1415926535897932
  8.  
  9. static const double PI = pi_,
  10.     twopi =pi_ * 2.0,
  11.     piover4= pi_ /    4.0    ,
  12.     fouroverpi= 4.0 / pi_ ,
  13.     piover2= pi_ /    2.0;
  14.  
  15. static double atanc[] =    {
  16.     0.0,
  17.     0.4636476090008061165,
  18.     0.7853981633974483094,
  19.     0.98279372324732906714,
  20.     1.1071487177940905022,
  21.     1.1902899496825317322,
  22.     1.2490457723982544262,
  23.     1.2924966677897852673,
  24.     1.3258176636680324644
  25. };
  26.  
  27. /*  aint(x)      Return integer part of number.  Truncates towards 0     */
  28. #define aint(x)    ((double)((int)(x)))
  29. #define    fabs(x)     ((x < 0.0) ? -x : x)
  30.  
  31.  
  32. /*  sin(x)      Return sine, x in radians  */
  33.  
  34. double sin(double x)
  35. {
  36.     int     sign;
  37.     double     y, r, z;
  38.  
  39.     x = ( (sign= (x    < 0.0) ) ?    -x: x);
  40.  
  41.     if (x >    twopi)
  42.        x -=    (aint(x    / twopi) * twopi);
  43.  
  44.     if (x >    PI) {
  45.        x -=    PI;
  46.        sign    = !sign;
  47.     }
  48.  
  49.     if (x >    piover2)
  50.        x = PI - x;
  51.  
  52.     if (x <    piover4) {
  53.        y = x * fouroverpi;
  54.        z = y * y;
  55.        r = y * (((((((-0.202253129293E-13 *    z + 0.69481520350522E-11) * z -
  56.           0.17572474176170806E-8) *    z + 0.313361688917325348E-6) * z -
  57.           0.365762041821464001E-4) * z + 0.249039457019271628E-2) *    z -
  58.           0.0807455121882807815) * z + 0.785398163397448310);
  59.     } else {
  60.        y = (piover2    - x) * fouroverpi;
  61.        z = y * y;
  62.        r = ((((((-0.38577620372E-12    * z + 0.11500497024263E-9) * z -
  63.           0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z    -
  64.           0.325991886926687550E-3) * z + 0.0158543442438154109) * z    -
  65.           0.308425137534042452) * z    + 1.0;
  66.     }
  67.     return sign ? -r : r;
  68. }
  69.  
  70.  
  71. double cos(double x)
  72. {
  73.     x = (x < 0.0) ?    -x : x;
  74.     if (x >    twopi)              /* Do range reduction here to limit */
  75.        x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2    */
  76.     return sin(x + piover2);
  77. }
  78.  
  79. /*  tan(x)      Return tangent, x in radians,    by identity  */
  80.  
  81. double tan(double x)
  82. {
  83.     return sin(x) /    cos(x);
  84. }
  85.  
  86. /*  sqrt(x)      Return square    root.  Initial guess, then Newton-
  87.           Raphson refinement  */
  88.  
  89. double sqrt(double x)
  90. {
  91.     double c, cl, y;
  92.     int n;
  93.  
  94.     if (x == 0.0)
  95.        return 0.0;
  96.  
  97.  
  98.     y = (0.154116 +    1.893872 * x) /    (1.0 + 1.047988    * x);
  99.  
  100.     c = (y - x / y)    / 2.0;
  101.     cl = 0.0;
  102.     for (n = 50; c != cl &&    n--;) {
  103.        y = y - c;
  104.        cl =    c;
  105.        c = (y - x /    y) / 2.0;
  106.     }
  107.     return y;
  108. }
  109.  
  110. /*  atan(x)      Return arctangent in radians,
  111.           range    -PI/2 to PI/2  */
  112.  
  113. double atan(double x)
  114. {
  115.     int sign, l, y;
  116.     double a, b, z;
  117.  
  118.     x = (((sign = (x < 0.0)) != 0) ? -x : x);
  119.     l = 0;
  120.  
  121.     if (x >= 4.0) {
  122.        l = -1;
  123.        x = 1.0 / x;
  124.        y = 0;
  125.        goto    atl;
  126.     } else {
  127.        if (x < 0.25) {
  128.           y    = 0;
  129.           goto atl;
  130.        }
  131.     }
  132.  
  133.     y = aint(x / 0.5);
  134.     z = y *    0.5;
  135.     x = (x - z) / (x * z + 1);
  136.  
  137. atl:
  138.     z = x *    x;
  139.     b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) *    z +
  140.         1277025750.0) * z +    1550674125.0) *    z + 654729075.0;
  141.     a = (((13852575.0 * z +    216602100.0) * z + 891080190.0)    * z +
  142.         1332431100.0) * z +    654729075.0;
  143.     a = (a / b) * x    + atanc[y];
  144.     if (l)
  145.        a=piover2 - a;
  146.     return sign ? -a : a;
  147. }
  148.  
  149. /*  atan2(y,x)      Return arctangent in radians of y/x,
  150.           range    -PI to PI  */
  151.  
  152. double atan2(double y, double x)
  153. {
  154.     double temp;
  155.  
  156.     if (x == 0.0) {
  157.        if (y == 0.0)   /*  Special case: atan2(0,0)    = 0  */
  158.           return 0.0;
  159.        else    if (y >    0)
  160.           return piover2;
  161.        else
  162.           return -piover2;
  163.     }
  164.     temp = atan(y /    x);
  165.     if (x <    0.0) {
  166.        if (y >= 0.0)
  167.           temp += pi_;
  168.        else
  169.           temp -= pi_;
  170.     }
  171.     return temp;
  172. }
  173.  
  174. /*  asin(x)      Return arcsine in radians of x  */
  175.  
  176. double asin(double x)
  177. {
  178.     return atan2(x,    sqrt(1 - x * x));
  179. }
  180.